home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / basic / ubmalm.zip / dldiosys.ub < prev    next >
Text File  |  1990-08-22  |  6KB  |  146 lines
  1.    10   ' Driver program for LDIOSYS
  2.    12   ' 19 June 1990
  3.    15   word 20
  4.    20   dim A(15,16),V(15,15)
  5.    30   print "How many rows in the coeff. matrix":input R
  6.    40   print "How many columns":input C
  7.    80   gosub *Getarray(R,C,&A())
  8.   100   gosub *Ldiosys(R,C,&A(),&V(),&Soln,&Rank)
  9.   290   print =print + lprint
  10.   300   gosub *Report(R,C,&A(),&V(),&Soln,&Rank)
  11.   302   print:print:print:print =lprint +"Malm.UBD"
  12.   305   gosub *Report2(R,C,&A(),&V(),&Soln,&Rank)
  13.   306   print:print:print:print:print:
  14.   310   print =print
  15.   400   end
  16.   490   '
  17.   495   '
  18.   500   *Getarray(R,C,&A())
  19.   510   ' Obtains from the keyboard the coefficient array A.
  20.   520   '15 June 1990
  21.   530   local I%,J%,T%
  22.   540   print "Enter the equation coefficients including the R.H.S."
  23.   550   print "Follow each number with return."
  24.   560   for I%=1 to R:for J%=1 to C+1:input A(I%,J%):next J%
  25.   570   print "eq'n completed":next I%
  26.   580   print "The matrix (by columns) is:"
  27.   590   for J%=1 to C+1:for I%=1 to R:print A(I%,J%):next I%
  28.   600   print "Lightly press return to continue":T%=input$(1):print
  29.   610   :next J%
  30.   620   repeat input "Do you want to change an entry (1 or 0)";T%
  31.   630   if T%=1 then
  32.   640   :print "Which row ":input I%:print "Which column ":input J%
  33.   650   :print "New value: ":input A(I%,J%)
  34.   660   :endif
  35.   670   until T%=0
  36.   680   return ' End of subroutine Getarray.
  37.   790   '
  38.   795   '
  39.   800   *Report(R,C,&A(),&V(),&Soln,&Rank)
  40.   810   ' Outputs the results in a format that is O.K. if the matrices are not
  41.   820   ' too big nor are the numbers.
  42.   830   ' 18 June 1990
  43.   850   local I%,J%,S
  44.   860   if Soln=0 then print "No solution." endif
  45.   870   if Rank=0 then print "The rank is zero." else
  46.   875   :print "The rank is ";Rank
  47.   880   :print "The augmented coeff. matrix followed by the record matrix:"
  48.   890   :for I%=1 to R:print "        ";
  49.   900   :for J%=1 to C+1:print A(I%,J%),
  50.   910   :next J%:print:next I%
  51.   920   :for I%=1 to C:print "        ";
  52.   930   :for J%=1 to C:print V(I%,J%),
  53.   940   :next J%:print:next I%
  54.   950   :for I%=1 to 3:print:next I%
  55.   960   :if Soln then
  56.   970   ::for I%=1 to C:S=0:print "        ";
  57.   980   ::for J%=1 to Rank:S+=V(I%,J%)*A(J%,C+1):next J%
  58.   990   ::print S,:print "          ";
  59.  1000   ::for J%=Rank+1 to C:print V(I%,J%),:next J%
  60.  1010   ::print:next I%
  61.  1020   ::endif ' if soln
  62.  1030   :endif:return ' End of SubroutineReport.
  63.  1090   '
  64.  1095   '
  65.  1100   *Report2(R,C,&A(),&V(),&Soln,&Rank)
  66.  1110   ' Outputs the results one number to a line.
  67.  1130   ' 19 June 1990.
  68.  1140   local I%,J%,S
  69.  1160   print "The augmented coefficient matrix by columns:":print
  70.  1170   for J%=1 to C+1:for I%=1 to R:print A(I%,J%):next I%:print:next J%
  71.  1200   print:print "The record matrix by columns:":print
  72.  1210   for J%=1 to C:for I%=1 to C:print V(I%,J%):next I%:print:next J%
  73.  1230   print:print "The rank is ";Rank
  74.  1240   if Soln=0 then print "No solution":return endif
  75.  1250   print "A particular solution (column):":print
  76.  1260   for I%=1 to C:S=0
  77.  1270   for J%=1 to Rank:S+=V(I%,J%)*A(J%,C+1):next J%
  78.  1280   print S:next I%:print
  79.  1300   print "The coeff. matrix of the parameters(by columns):":print
  80.  1310   for I%=Rank+1 to C:for J%=1 to C:print V(J%,I%):next J%:print:next I%
  81.  1320   return ' End of subroutine Report2.
  82.  1330   '
  83.  1335   '
  84.  2320   *LDiosys(R,C,&A(),&V(),&Soln,&Rank)
  85.  2330   ' Solve the linear Diophantine system represented by the (augmented)
  86.  2340   ' matrix A().  Matrix V() contains a record of the manipulations.
  87.  2350   ' A is R by C+1, while V is C by C.  Soln is 1 if there is a solution
  88.  2360   ' and is 0 if there is no solution.  Rank is the rank of the system.
  89.  2370   '
  90.  2380   ' NOTE - there is a GLOBAL variable J% - required to receive the value
  91.  2390   ' of the function Pivotind, which LDiosys calls.
  92.  2400   '
  93.  2410   ' 19 June 1990
  94.  2420   local K%,M%,N%,I%,Done,Te,Found,P
  95.  2430   for I%=1 to C:for K%=1 to C:V(I%,K%)=0:next K%:next I%
  96.  2440   for I%=1 to C:V(I%,I%)=1:next I%
  97.  2450   Done=0:I%=0:Soln=1
  98.  2460   while and{I%<R,I%<C,Done=0}
  99.  2470   inc I%:J%=fnPivotind(I%)
  100.  2480   if J%=0 then K%=I%:Found=0
  101.  2490   :while and{Found=0,K%<R}
  102.  2500   :inc K%:J%=fnPivotind(K%)
  103.  2510   :if J%<>0 then Found=1 endif
  104.  2520   :wend
  105.  2530   :if Found=0 then Done=1 else
  106.  2540   :for M%=1 to C+1:swap A(K%,M%),A(I%,M):next M%
  107.  2550   :endif endif
  108.  2560   if Done=0 then
  109.  2570   :repeat K%=J%:P=A(I%,J%)
  110.  2580   :for N%=I% to C
  111.  2590   :if N%<>J% then Te=A(I%,N%)\P
  112.  2600   :for M%=I% to R:A(M%,N%)-=Te*A(M%,J%):next M%
  113.  2610   :for M%=1 to C:V(M%,N%)-=Te*V(M%,J%):next M%
  114.  2620   :endif:next N%
  115.  2630   :J%=fnPivotind(I%)
  116.  2640   :until J%=K%
  117.  2650   :if I%<>J% then
  118.  2660   :for M%=I% to R:swap A(M%,I%),A(M%,J%):next M%
  119.  2670   :for M%=1 to C:swap V(M%,I%),V(M%,J%):next M%
  120.  2680   :endif
  121.  2690   :if A(I%,C+1)@A(I%,I%)<>0 then Done=1 else
  122.  2700   :A(I%,C+1)\=A(I%,I%):A(I%,I%)=1
  123.  2710   :for M%=I%+1 to R:A(M%,C+1)-=A(M%,I%)*A(I%,C+1):A(M%,I%)=0:next M%
  124.  2720   :endif
  125.  2730   :endif ' done=0
  126.  2740   wend
  127.  2750   Te=min(R,C):Rank=0
  128.  2760   for K%=1 to Te:if A(K%,K%)<>0 then inc Rank endif next K%
  129.  2770   K%=0
  130.  2780   while and{K%<Rank,Soln}:inc K%:if A(K%,K%)<>1 then Soln=0 endif wend
  131.  2790   for K%=Rank+1 to R:if A(K%,C+1)<>0 then Soln=0 endif:next K%
  132.  2800   return ' End of subroutine LDiosys.
  133.  2810   '
  134.  2820   fnPivotind(Kk%)
  135.  2830   ' Returns the column index of the smallest (in absolute value )
  136.  2840   ' non-zero entry in row Kk% of the coefficient matrix A.
  137.  2850   ' If all entries are zero, then Pivotind = 0.
  138.  2855   ' 15 June 1990
  139.  2860   local T,S,K%,J%
  140.  2870   T=abs(A(Kk%,1))
  141.  2880   if T<>0 then K%=1 else K%=0 endif
  142.  2890   for J%=2 to C
  143.  2900   S=abs(A(Kk%,J%))
  144.  2910   if and{S<>0,or{T=0,S<T}} then T=S:K%=J% endif
  145.  2920   next J%:return(K%) ' End of function Pivotind.
  146.